germination
library(rethinking)Loading required package: rstan
Loading required package: StanHeaders
Loading required package: ggplot2
Keep up to date with changes at https://www.tidyverse.org/blog/
rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Loading required package: parallel
rethinking (Version 2.11)
Attaching package: ‘rethinking’
The following object is masked from ‘package:stats’:
rstudentlibrary(brms)Loading required package: Rcpp
Loading 'brms' package (version 2.13.5). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: ‘brms’
The following objects are masked from ‘package:rethinking’:
LOO, stancode, WAIC
The following object is masked from ‘package:rstan’:
loo
The following object is masked from ‘package:stats’:
arlibrary(tidyverse)Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble 3.0.3 ✓ dplyr 1.0.0
✓ tidyr 1.1.0 ✓ stringr 1.4.0
✓ readr 1.3.1 ✓ forcats 0.5.0
✓ purrr 0.3.4
── Conflicts ───────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::extract() masks rstan::extract()
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
x purrr::map() masks rethinking::map()germ <- read_csv("../Dimensions/hydrothermal-all-species/data/light_round1_tall.csv") %>%
filter(wps == 0) %>%
select(pops, temps, total_seeds, germ, day, cumulative_germ)Parsed with column specification:
cols(
pops = col_character(),
temps = col_double(),
wps = col_double(),
date = col_character(),
total_seeds = col_double(),
germ = col_double(),
start_date = col_date(format = ""),
census_date = col_date(format = ""),
day = col_double(),
cumulative_germ = col_double(),
cumulative_prop_germ = col_double()
)germwhat if need one event per row:
one_per_row <- function(df) {
total_seed <- max(df$total_seeds, sum(df$germ))
newdata <- tibble(id=1:total_seed, germ=0, day=max(df$day))
df <- df %>% filter(germ>0)
count <- 1
if (nrow(df) > 0) {
for (i in 1:nrow(df)) { # we look at each row of the df where germination occured
for (j in 1:df$germ[i]) { # now update the newdata to reflect the germiantion of each seed
newdata$germ[count] <- 1
newdata$day[count]=df$day[i]
count <- count+1 # count keeps track of which individual we are at in the new data
} # for j
} # for i
} # if
return(newdata)
}
germone <- germ %>% group_by(pops, temps) %>%
select(-cumulative_germ) %>% # not needed in this encoding (I think...in any case would need to be recalculated)
nest() %>%
mutate(newdata=map(data, one_per_row)) %>%
select(-data) %>%
unnest(newdata)
germone1
Is effect of temp ~linear on cum germ?
germ %>% filter(pops=="STDI", day==28) %>%
ggplot(aes(x=temps,y=cumulative_germ)) +
geom_col()
no, so we can’t treat temp as a continuous linear predictor
germ.stdi <- germone %>% filter(pops=="STDI") %>% dplyr::select(-pops)Adding missing grouping variables: `pops`germ.stdi2
censored exponential
rethinking
exponential rate curve, censoring seed that don’t germinate.
d <- list(germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
m1.1 <- ulam(
alist(
day | germ==1 ~ exponential( lambda),
day | germ==0 ~ custom(exponential_lccdf( !Y | lambda)),
lambda <- 1.0 / mu,
log(mu) <- a[temps],
a[temps] ~ normal(0,1)),
data=d,
chains=4,
cores = 4
)starting worker pid=6368 on localhost:11907 at 09:03:25.275
starting worker pid=6382 on localhost:11907 at 09:03:25.610
starting worker pid=6397 on localhost:11907 at 09:03:25.887
starting worker pid=6412 on localhost:11907 at 09:03:26.197
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 8.2e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 8.6e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.86 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.223272 seconds (Warm-up)
Chain 1: 0.22504 seconds (Sampling)
Chain 1: 0.448312 seconds (Total)
Chain 1:
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 7.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.278519 seconds (Warm-up)
Chain 2: 0.208053 seconds (Sampling)
Chain 2: 0.486572 seconds (Total)
Chain 2:
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.274003 seconds (Warm-up)
Chain 3: 0.210096 seconds (Sampling)
Chain 3: 0.484099 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8.4e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.84 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.260888 seconds (Warm-up)
Chain 4: 0.217804 seconds (Sampling)
Chain 4: 0.478692 seconds (Total)
Chain 4: precis(m1.1, depth = 2)The above represent log(mean time to germination)
exp(5.45)[1] 232.7582exp(2.48)[1] 11.94126posterior
preddata <- expand_grid(temps=1:8, day=1:28)
pred <- link(m1.1, data = preddata)
str(pred)List of 2
$ lambda: num [1:2000, 1:224] 0.00585 0.00443 0.00493 0.0032 0.00844 ...
$ mu : num [1:2000, 1:224] 171 226 203 312 119 ...mu is average days to germ, lambda is rate parameter. neither change over time, of course, so having days doesn’t really make sense.
preddata$mu <- apply(pred$mu,2,mean)
preddata$low <- apply(pred$mu,2,HPDI)[1,]
preddata$high <- apply(pred$mu, 2, HPDI)[2,]Single temp. Don’t need day posterior
preddata <- expand_grid(temps=3)
pred <- link(m1.1, data = preddata)
str(pred)how to convert to probs? use pexp.
predprobs <- pexp(1:28,rate=pred$lambda[1])#plot(x=1:28,y=predprobs, type="l") # crashes on ploteven though it isn’t using day, including day in the prediction data frame will help me keep data in the correct format. Maybe.
preddata <- expand_grid(temps=1:8, day=1:28)
pred <- link(m1.1, data = preddata)
str(pred)List of 2
$ lambda: num [1:2000, 1:224] 0.00585 0.00443 0.00493 0.0032 0.00844 ...
$ mu : num [1:2000, 1:224] 171 226 203 312 119 ...predresults <- preddata %>%
mutate(lambda=as.list(as.data.frame(pred$lambda)))
predresultspredresults <- predresults %>%
mutate(probs=map2(day, lambda, ~ pexp(.x, .y)),
mu=map_dbl(probs, mean),
lower=map_dbl(probs, ~ HPDI(.)[1] %>% unlist()),
upper=map_dbl(probs, ~ HPDI(.)[2]) %>% unlist())
predresultspredresults %>% select(-lambda, -probs) %>%
mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
geom_line() 
Add realdata:
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
predresults %>% select(-lambda, -probs) %>%
mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
geom_line() +
geom_point(aes(y=prop_germ), data=stdi.plot)
Poor!
brms censoring model
need to set up indicator for censoring.
germ.stdi <- germ.stdi %>%
mutate(cens=ifelse(germ==0, "right", "none"),
tempsc=as.character(temps) %>% str_pad(width=2, pad="0"))
germ.stdiget_prior(day | cens(cens) ~ 0 + tempsc, family = exponential, data=germ.stdi)m1.2 <- brm(day | cens(cens) ~ 0 + tempsc,
family = exponential(),
set_prior("normal(0,1)", class="b"),
data = germ.stdi, )Compiling Stan program...
Start sampling
SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 9.8e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.423825 seconds (Warm-up)
Chain 1: 0.429479 seconds (Sampling)
Chain 1: 0.853304 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 6.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.426314 seconds (Warm-up)
Chain 2: 0.455533 seconds (Sampling)
Chain 2: 0.881847 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 6.2e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.43039 seconds (Warm-up)
Chain 3: 0.41959 seconds (Sampling)
Chain 3: 0.84998 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000196 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.96 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.430918 seconds (Warm-up)
Chain 4: 0.426936 seconds (Sampling)
Chain 4: 0.857854 seconds (Total)
Chain 4: summary(m1.2) Family: exponential
Links: mu = log
Formula: day | cens(cens) ~ 0 + tempsc
Data: germ.stdi (Number of observations: 396)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tempsc05 5.44 0.37 4.76 6.22 1.00 6689 2848
tempsc10 4.70 0.28 4.18 5.30 1.00 7327 2780
tempsc15 4.02 0.23 3.60 4.49 1.00 7811 2805
tempsc20 2.48 0.16 2.17 2.81 1.00 7247 2911
tempsc25 3.07 0.17 2.74 3.43 1.00 7072 3064
tempsc30 3.73 0.20 3.36 4.14 1.00 7913 3195
tempsc35 5.21 0.35 4.57 5.94 1.00 6135 3096
tempsc40 5.62 0.40 4.89 6.46 1.00 7016 2782
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).predict(m1.2) %>% head Estimate Est.Error Q2.5 Q97.5
[1,] 249.2430 286.1171 6.667386 1031.5988
[2,] 242.8674 286.0824 5.156902 996.2780
[3,] 252.4897 300.5946 5.439334 1066.8564
[4,] 247.2513 284.7164 5.992272 1018.9538
[5,] 251.3389 286.2790 5.726633 1042.6966
[6,] 243.3522 263.6542 6.032181 953.4948log(mean time to germination)
predict(m1.2, newdata = expand_grid(tempsc=unique(germ.stdi$tempsc), cens=c("none", "right"))) Estimate Est.Error Q2.5 Q97.5
[1,] 250.54096 285.42410 5.4041250 1048.05543
[2,] 249.55624 296.61765 5.9321050 1064.49475
[3,] 115.14309 125.64755 2.5715721 458.46359
[4,] 114.02507 119.67649 2.7638015 441.18807
[5,] 57.28506 59.16199 1.4965791 216.83994
[6,] 57.09426 58.62516 1.4054719 218.42622
[7,] 12.12453 12.28691 0.3618226 47.03079
[8,] 12.19735 12.56483 0.3080198 45.58575
[9,] 21.67803 22.76933 0.5948858 82.15578
[10,] 21.86662 21.88403 0.5093664 78.33396
[11,] 42.78315 43.77602 1.0944283 161.49532
[12,] 43.82032 45.09836 1.1251890 162.88073
[13,] 192.67635 215.60433 3.8413475 773.76129
[14,] 189.78430 217.32587 4.2042633 784.30106
[15,] 300.67681 347.99351 6.5632012 1292.12783
[16,] 293.24047 343.31921 7.0867706 1252.76443precis(m1.1, depth = 2)again, these are log(mean time to germination)
lambda should be 1/mu, so
#plot(1:28,pexp(1:28, 1/exp(2.48)), type="l", col="red")
#lines(1:28,pexp(1:28, 1/exp(5.60)), type="l", col="blue")
#crashesget_prior( day | cens(cens) ~ 0 + tempsc,
family = Gamma(),
data = germ.stdi)m1.3 <- brm(day | cens(cens) ~ 0 + tempsc,
family = Gamma(),
set_prior("normal(0,1)", class="b", lb=0),
data = germ.stdi)Compiling Stan program...
Start sampling
SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000441 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.41 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 3.83018 seconds (Warm-up)
Chain 1: 2.69174 seconds (Sampling)
Chain 1: 6.52191 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000423 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.23 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 3.7561 seconds (Warm-up)
Chain 2: 2.59898 seconds (Sampling)
Chain 2: 6.35508 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000369 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.69 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 3.65055 seconds (Warm-up)
Chain 3: 2.61099 seconds (Sampling)
Chain 3: 6.26154 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000518 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.18 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 3.6174 seconds (Warm-up)
Chain 4: 2.67537 seconds (Sampling)
Chain 4: 6.29276 seconds (Total)
Chain 4: summary(m1.3) Family: gamma
Links: mu = inverse; shape = identity
Formula: day | cens(cens) ~ 0 + tempsc
Data: germ.stdi (Number of observations: 396)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tempsc05 0.00 0.00 0.00 0.00 1.00 5074 2617
tempsc10 0.00 0.00 0.00 0.01 1.00 4408 2321
tempsc15 0.01 0.00 0.00 0.02 1.00 4649 2755
tempsc20 0.07 0.02 0.05 0.11 1.00 5340 2742
tempsc25 0.04 0.01 0.02 0.06 1.00 5426 2802
tempsc30 0.02 0.00 0.01 0.03 1.00 5425 2945
tempsc35 0.00 0.00 0.00 0.00 1.00 4444 2688
tempsc40 0.00 0.00 0.00 0.00 1.00 4553 2247
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.68 0.07 0.56 0.82 1.00 3100 2848
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).newdata <- expand_grid(tempsc=unique(germ.stdi$tempsc), cens=c("none", "right"))
cbind(newdata, predict(m1.3, newdata = newdata))The problem with these models is that they assume that the censored observations are controlled by the same process that is controlling the ones that do germinate. I do not think that is correct, we need a mixture model.
Trying a mixture model for a ZI exponential
look at ulam censored exponential
rethinking::stancode(m1.1)data{
int germ[396];
vector[396] day;
int temps[396];
}
parameters{
vector[8] a;
}
model{
vector[396] lambda;
vector[396] mu;
a ~ normal( 0 , 1 );
for ( i in 1:396 ) {
mu[i] = a[temps[i]];
mu[i] = exp(mu[i]);
}
for ( i in 1:396 ) {
lambda[i] = 1/mu[i];
}
for ( i in 1:396 )
if ( germ[i] == 0 ) target += exponential_lccdf(day[i] | lambda[i]);
for ( i in 1:396 )
if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `day`. look at brms censored exponential
stancode(m1.2)// generated with brms 2.13.5
functions {
}
data {
int<lower=1> N; // number of observations
vector[N] Y; // response variable
int<lower=-1,upper=2> cens[N]; // indicates censoring
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
}
parameters {
vector[K] b; // population-level effects
}
transformed parameters {
}
model {
// initialize linear predictor term
vector[N] mu = X * b;
for (n in 1:N) {
// apply the inverse link function
mu[n] = exp(-(mu[n]));
}
// priors including all constants
target += normal_lpdf(b | 0,1);
// likelihood including all constants
if (!prior_only) {
for (n in 1:N) {
// special treatment of censored data
if (cens[n] == 0) {
target += exponential_lpdf(Y[n] | mu[n]);
} else if (cens[n] == 1) {
target += exponential_lccdf(Y[n] | mu[n]);
} else if (cens[n] == -1) {
target += exponential_lcdf(Y[n] | mu[n]);
}
}
}
}
generated quantities {
}Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `day`.
3: Unknown or uninitialised column: `germ`.
4: Unknown or uninitialised column: `day`. Zero inflated, as generated by ulam:
There were 18 warnings (use warnings() to see them)# define parameters
prob_drink <- 0.2 # 20% of days
rate_work <- 1 # average 1 manuscript per day
# sample one year of production
N <- 365
# simulate days monks drink
set.seed(365)
drink <- rbinom( N , 1 , prob_drink )
# simulate manuscripts completed
y <- (1-drink)*rpois( N , rate_work )
rethinking::stancode(ulam( alist(
y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
logit(p) <- ap,
log(lambda) <- al,
ap ~ dnorm(-1.5,1),
al ~ dnorm(1,0.5)),
data=list(y=as.integer(y)) , sample=FALSE))data{
int y[365];
}
parameters{
real ap;
real al;
}
model{
real p;
real lambda;
al ~ normal( 1 , 0.5 );
ap ~ normal( -1.5 , 1 );
lambda = al;
lambda = exp(lambda);
p = ap;
p = inv_logit(p);
for ( i in 1:365 )
if ( y[i] == 0 ) target += log_mix(p, 0, poisson_lpmf(0 | lambda));
for ( i in 1:365 )
if ( y[i] > 0 ) target += log1m(p) + poisson_lpmf(y[i] | lambda);
}Zero inflated exponential model, attempt 1
d <- list(N=nrow(germ.stdi),Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `day`.
3: Unknown or uninitialised column: `germ`.
4: Unknown or uninitialised column: `day`.
5: Unknown or uninitialised column: `germ`.
6: Unknown or uninitialised column: `day`. germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.4 <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
vector[8] a; // alphas for the linear model
real ap; // probability of drink
}
model{
real p;
vector[N] lambda;
vector[N] mu;
a ~ normal( 0 , 1 );
ap ~ normal( 0, 1 );
p = ap;
p = inv_logit(p); // inverse link function for probability drink
for ( i in 1:N ) {
mu[i] = a[temps[i]]; // linear model for mu
mu[i] = exp(mu[i]); // inverse link function part 1
lambda[i] = 1/mu[i]; // inverse link function part 2
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p, 0, exponential_lccdf(day[i] | lambda[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p) + exponential_lpdf(day[i] | lambda[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.4 <- stan(model_code=stanmodel1.4, data=d)
SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000112 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.12 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.769353 seconds (Warm-up)
Chain 1: 0.757463 seconds (Sampling)
Chain 1: 1.52682 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000103 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.855376 seconds (Warm-up)
Chain 2: 0.729975 seconds (Sampling)
Chain 2: 1.58535 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 9.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.95 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.798865 seconds (Warm-up)
Chain 3: 0.741906 seconds (Sampling)
Chain 3: 1.54077 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 9.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.801264 seconds (Warm-up)
Chain 4: 0.71224 seconds (Sampling)
Chain 4: 1.5135 seconds (Total)
Chain 4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-essTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-essprecis(m1.4, depth = 2)There were 16 warnings (use warnings() to see them)Model 1.4 fits a single dormancy rate for all, but probably we need to have separate dormancy rate for each temp
take a look at stancode for a binomial model to see what this looks like
db <- list(cumulative_germ=germ$cumulative_germ, temps=as.numeric(as.factor(germ$temps)))There were 12 warnings (use warnings() to see them)rethinking::stancode(ulam(alist(
cumulative_germ ~ dbern(mu),
mu <- logit(a + b[temps]),
a ~ dnorm(0,1),
b[temps] ~ dnorm(0,1)),
data = db, sample=FALSE, log_lik = TRUE))data{
int cumulative_germ[3840];
int temps[3840];
}
parameters{
real a;
vector[8] b;
}
model{
vector[3840] mu;
b ~ normal( 0 , 1 );
a ~ normal( 0 , 1 );
for ( i in 1:3840 ) {
mu[i] = logit(a + b[temps[i]]);
}
cumulative_germ ~ bernoulli( mu );
}
generated quantities{
vector[3840] log_lik;
vector[3840] mu;
for ( i in 1:3840 ) {
mu[i] = logit(a + b[temps[i]]);
}
for ( i in 1:3840 ) log_lik[i] = bernoulli_lpmf( cumulative_germ[i] | mu[i] );
}There is no inverse logit because it is built in to the bernoulli function
Try to write a model with different dormancy probs for the different temps.
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.5 <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
vector[8] a; // alpha for the exponential curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] lambda;
vector[N] mu;
a ~ normal( 0 , 1 );
ap ~ normal( 0, 1 );
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]); // inverse link function
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
mu[i] = exp(mu[i]); // inverse link function
}
for ( i in 1:N ) {
lambda[i] = 1/mu[i];
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, exponential_lccdf(day[i] | lambda[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + exponential_lpdf(day[i] | lambda[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.5 <- stan(model_code=stanmodel1.5, data=d)
SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000202 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1.21524 seconds (Warm-up)
Chain 1: 0.957846 seconds (Sampling)
Chain 1: 2.17309 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000121 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.21 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.16162 seconds (Warm-up)
Chain 2: 0.969372 seconds (Sampling)
Chain 2: 2.13099 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000152 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.18377 seconds (Warm-up)
Chain 3: 0.939891 seconds (Sampling)
Chain 3: 2.12366 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.00012 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 1.17349 seconds (Warm-up)
Chain 4: 0.859005 seconds (Sampling)
Chain 4: 2.03249 seconds (Total)
Chain 4: precis(m1.5, depth = 2)shinystan::launch_shinystan(m1.5)
Launching ShinyStan interface... for large models this may take some time.
Listening on http://127.0.0.1:3321post <- as.data.frame(m1.5)
head(post)post.cor <- cor(post) %>% round(2)
post.cor[abs(post.cor)<.25] <- NA
post.cor a[1] a[2] a[3] a[4] a[5] a[6] a[7] a[8] ap[1] ap[2] ap[3] ap[4] ap[5] ap[6]
a[1] 1.00 NA NA NA NA NA NA NA -0.39 NA NA NA NA NA
a[2] NA 1 NA NA NA NA NA NA NA NA NA NA NA NA
a[3] NA NA 1 NA NA NA NA NA NA NA NA NA NA NA
a[4] NA NA NA 1 NA NA NA NA NA NA NA NA NA NA
a[5] NA NA NA NA 1 NA NA NA NA NA NA NA NA NA
a[6] NA NA NA NA NA 1 NA NA NA NA NA NA NA NA
a[7] NA NA NA NA NA NA 1.00 NA NA NA NA NA NA NA
a[8] NA NA NA NA NA NA NA 1 NA NA NA NA NA NA
ap[1] -0.39 NA NA NA NA NA NA NA 1.00 NA NA NA NA NA
ap[2] NA NA NA NA NA NA NA NA NA 1 NA NA NA NA
ap[3] NA NA NA NA NA NA NA NA NA NA 1 NA NA NA
ap[4] NA NA NA NA NA NA NA NA NA NA NA 1 NA NA
ap[5] NA NA NA NA NA NA NA NA NA NA NA NA 1 NA
ap[6] NA NA NA NA NA NA NA NA NA NA NA NA NA 1
ap[7] NA NA NA NA NA NA -0.32 NA NA NA NA NA NA NA
ap[8] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
lp__ NA NA NA NA NA NA NA NA NA NA NA NA NA NA
ap[7] ap[8] lp__
a[1] NA NA NA
a[2] NA NA NA
a[3] NA NA NA
a[4] NA NA NA
a[5] NA NA NA
a[6] NA NA NA
a[7] -0.32 NA NA
a[8] NA NA NA
ap[1] NA NA NA
ap[2] NA NA NA
ap[3] NA NA NA
ap[4] NA NA NA
ap[5] NA NA NA
ap[6] NA NA NA
ap[7] 1.00 NA NA
ap[8] NA 1 NA
lp__ NA NA 1now need to make predictions
post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu_lambda <- post_logavgdays %>%
summarise_all(mean) %>% # posterior mean
mutate_all(~ 1/exp(.)) # convert to lambda (inverse link)
mu_lambdapost_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_pposterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
lambda=t(mu_lambda),
p=t(mu_p))plot, include original data
post_plot <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pexp(day, rate=lambda) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAadjust priors
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.5a <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
vector[8] a; // alpha for the exponential curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] lambda;
vector[N] mu;
a ~ normal( 0 , 2 );
ap ~ normal( 0, 2 );
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]); // inverse link function
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
mu[i] = exp(mu[i]); // inverse link function
}
for ( i in 1:N ) {
lambda[i] = 1/mu[i];
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, exponential_lccdf(day[i] | lambda[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + exponential_lpdf(day[i] | lambda[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.5a <- stan(model_code=stanmodel1.5a, data=d)
SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00016 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.6 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 1.2908 seconds (Warm-up)
Chain 1: 0.945991 seconds (Sampling)
Chain 1: 2.23679 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000136 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.36 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 1.49546 seconds (Warm-up)
Chain 2: 1.02455 seconds (Sampling)
Chain 2: 2.52001 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000121 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.21 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.53961 seconds (Warm-up)
Chain 3: 1.05771 seconds (Sampling)
Chain 3: 2.59732 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.00012 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 1.43043 seconds (Warm-up)
Chain 4: 1.27791 seconds (Sampling)
Chain 4: 2.70834 seconds (Total)
Chain 4: There were 1 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problemsprecis(m1.5a, depth = 2)shinystan::launch_shinystan(m1.5a)post <- as.data.frame(m1.5a)
post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu_lambda <- post_logavgdays %>%
summarise_all(mean) %>% # posterior mean
mutate_all(~ 1/exp(.)) # convert to lambda (inverse link)
mu_lambda
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
lambda=t(mu_lambda),
p=t(mu_p))
post_plot <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pexp(day, rate=lambda) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAstill not great
gamma
Look at gamma stancode from brms:
make_stancode(day ~ tempsc,
family="Gamma",
prior=set_prior("normal(0,1)"),
data=germ.stdi,
sample=FALSE)// generated with brms 2.13.5
functions {
}
data {
int<lower=1> N; // number of observations
vector[N] Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
real<lower=0> shape; // shape parameter
}
transformed parameters {
}
model {
// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
for (n in 1:N) {
// apply the inverse link function
mu[n] = shape * exp(-(mu[n]));
}
// priors including all constants
target += normal_lpdf(b | 0,1);
target += student_t_lpdf(Intercept | 3, 3.3, 2.5);
target += gamma_lpdf(shape | 0.01, 0.01);
// likelihood including all constants
if (!prior_only) {
target += gamma_lpdf(Y | shape, mu);
}
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}Try developing a ZI gamma model:
stan uses that gamma(alpha, beta) parameterization where “alpha” is shape" and “beta” is “rate”. So try having the shape be the same for all temps.
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.7 <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
real<lower=0> shape; // should set lower bound
vector[8] a; // alpha for the gamma curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] mu;
a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
ap ~ normal( 0, 1.5 );
shape ~ normal(0, .5); //narrow priors to overcome divergent transitions
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]);
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
// apply the inverse link function
mu[i] = shape * exp(-(mu[i]));
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7 <- stan(model_code=stanmodel1.7, data=d, chains=4, cores=4, control=list(adapt_delta=.99))starting worker pid=95474 on localhost:11147 at 12:57:13.894
starting worker pid=95488 on localhost:11147 at 12:57:14.340
starting worker pid=95503 on localhost:11147 at 12:57:14.749
starting worker pid=95518 on localhost:11147 at 12:57:15.343
SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.001004 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.04 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000521 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.21 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.00068 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.8 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000521 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.21 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 11.1101 seconds (Warm-up)
Chain 1: 7.85728 seconds (Sampling)
Chain 1: 18.9674 seconds (Total)
Chain 1:
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 11.5728 seconds (Warm-up)
Chain 2: 7.93563 seconds (Sampling)
Chain 2: 19.5084 seconds (Total)
Chain 2:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 12.2292 seconds (Warm-up)
Chain 4: 8.7347 seconds (Sampling)
Chain 4: 20.9639 seconds (Total)
Chain 4:
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 11.6817 seconds (Warm-up)
Chain 3: 14.9809 seconds (Sampling)
Chain 3: 26.6626 seconds (Total)
Chain 3: precis(m1.7, depth = 2)shinystan::launch_shinystan(m1.7)post <- as.data.frame(m1.7)
head(post)now need to make predictions
post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_ratemu_shape <- mean(post$shape)
mu_shape[1] 2.350936post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_pplot predictions and original data
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))
post_plot1.7 <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot1.7 %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAplaying with priors:
Make shape exponential
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.7a <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
real<lower=0> shape; // should set lower bound
vector[8] a; // alpha for the gamma curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] mu;
a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
ap ~ normal( 0, 1.5 );
shape ~ exponential(.5);
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]);
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
// apply the inverse link function
mu[i] = shape * exp(-(mu[i]));
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7a <- stan(model_code=stanmodel1.7a, data=d, chains=4, cores=4, control=list(adapt_delta=.99))starting worker pid=8458 on localhost:11686 at 09:35:49.473
starting worker pid=8472 on localhost:11686 at 09:35:49.776
starting worker pid=8487 on localhost:11686 at 09:35:50.142
starting worker pid=8502 on localhost:11686 at 09:35:50.398
SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Gradient evaluated at the initial value is not finite.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.000855 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000976 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.76 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000618 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000481 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.81 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 10.5168 seconds (Warm-up)
Chain 1: 7.28403 seconds (Sampling)
Chain 1: 17.8008 seconds (Total)
Chain 1:
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 11.3064 seconds (Warm-up)
Chain 3: 7.1714 seconds (Sampling)
Chain 3: 18.4778 seconds (Total)
Chain 3:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 14.172 seconds (Warm-up)
Chain 4: 7.23646 seconds (Sampling)
Chain 4: 21.4084 seconds (Total)
Chain 4:
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 10.5926 seconds (Warm-up)
Chain 2: 14.0434 seconds (Sampling)
Chain 2: 24.636 seconds (Total)
Chain 2: precis(m1.7a, depth = 2)post <- as.data.frame(m1.7a)
post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_ratemu_shape <- mean(post$shape)
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))
post_plot1.7a <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot1.7a %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAShape parameter is larger than in 1.7. The fit is slightly better at least for 20 and 25.
make a comparison plot
post_plot_compare <- bind_rows(old=post_plot1.7, new=post_plot1.7a, .id = "fit")
post_plot_compare %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=str_c(temps,fit),linetype=fit)) +
geom_line() +
geom_point(data=stdi.plot, aes(x=day,y=prop_germ,color=temps), inherit.aes = FALSE)
NAJust a hair better…
continue to work with priors
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.7b <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
real<lower=0> shape; // should set lower bound
vector[8] a; // alpha for the gamma curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] mu;
a ~ normal( 0 , .75 ); // divergent with 1
ap ~ normal( 0, 1.5 );
shape ~ exponential(1.5); // divergent if this is .5 and a is 0, 1
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]);
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
// apply the inverse link function
mu[i] = shape * exp(-(mu[i]));
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7b <- stan(model_code=stanmodel1.7b, data=d, chains=4, cores=4, control=list(adapt_delta=.99))starting worker pid=8771 on localhost:11686 at 09:42:59.153
starting worker pid=8785 on localhost:11686 at 09:42:59.459
starting worker pid=8800 on localhost:11686 at 09:42:59.791
starting worker pid=8815 on localhost:11686 at 09:43:00.136
SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000634 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.34 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000827 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 8.27 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000618 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.00093 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 9.3 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 12.1323 seconds (Warm-up)
Chain 1: 7.61628 seconds (Sampling)
Chain 1: 19.7486 seconds (Total)
Chain 1:
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 15.0125 seconds (Warm-up)
Chain 4: 7.62409 seconds (Sampling)
Chain 4: 22.6365 seconds (Total)
Chain 4:
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 12.2946 seconds (Warm-up)
Chain 3: 12.3514 seconds (Sampling)
Chain 3: 24.6461 seconds (Total)
Chain 3:
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 12.6255 seconds (Warm-up)
Chain 2: 14.3954 seconds (Sampling)
Chain 2: 27.0209 seconds (Total)
Chain 2: There were 39 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problemsprecis(m1.7b, depth = 2)post <- as.data.frame(m1.7b)
post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_ratemu_shape <- mean(post$shape)
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))
post_plot <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAno improvement over 1.7a
brms hurdle gamma
get_prior(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
family = hurdle_gamma,
data=germ.stdi)make_stancode(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
family = hurdle_gamma,
prior=set_prior("normal(0,1)"),
data=germ.stdi)m1.8 <- brm(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
family = hurdle_gamma,
prior=c(set_prior("normal(0,1)"),
set_prior("normal(0,1)", dpar="hu")), # must specify this!
data=germ.stdi, cores = 4)m1.8post <- as.data.frame(m1.8)
head(post)weird fit!
now need to make predictions
post_logavgdays <- post %>% select(shape, starts_with("b_temp"
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)post_logitp <- post %>% select(starts_with("b_hu")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_pposterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))plot predictions and original data
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
allow both shape and rate of ZI gamma to be predicted by temp
Don’t run. Doesn’t sample well because shape and rate are correlated.
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.9 <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
vector[8] as; // alpha for gama shape
vector[8] ar; // alpha for the gamma rate, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
vector[N] p;
vector[N] shape;
vector[N] rate;
as ~ normal( 0, 1 );
ar ~ normal( 0, 1 );
ap ~ normal( 0, 1.5 );
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]);
}
for ( i in 1:N ) {
shape[i] = as[temps[i]];
rate[i] = ar[temps[i]];
// apply the inverse link functions
shape[i] = exp(shape[i]);
rate[i] = exp(rate[i]);
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape[i], rate[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape[i], rate[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ gamma( shape[i], rate [i] );
}
generated quantities {
}
"
m1.9 <- stan(model_code=stanmodel1.9, data=d, chains=4, cores=4, control=list(adapt_delta=.99))precis(m1.9, depth = 2)post <- as.data.frame(m1.9)
post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))
post_plot <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
log normal
make_stancode(day ~ tempsc, family = lognormal(), data=germ.stdi)
// generated with brms 2.13.5
functions {
}
data {
int<lower=1> N; // number of observations
vector[N] Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
real<lower=0> sigma; // residual SD
}
transformed parameters {
}
model {
// initialize linear predictor term
vector[N] mu = Intercept + Xc * b;
// priors including all constants
target += student_t_lpdf(Intercept | 3, 3.3, 2.5);
target += student_t_lpdf(sigma | 3, 0, 2.5)
- 1 * student_t_lccdf(0 | 3, 0, 2.5);
// likelihood including all constants
if (!prior_only) {
target += lognormal_lpdf(Y | mu, sigma);
}
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.10 <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
vector[8] a; // alpha for the lognormal curve, one for each temp
vector[8] ap; // alpha for the proportion dormant, one for each temp
real<lower=0> sigma;
}
model{
vector[N] p;
vector[N] mu;
a ~ normal( 0 , 2 );
sigma ~ exponential (1);
ap ~ normal( 0, 2 );
for (i in 1:N) {
p[i] = ap[temps[i]];
p[i] = inv_logit(p[i]); // inverse link function
}
for ( i in 1:N ) {
mu[i] = a[temps[i]];
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, lognormal_lccdf(day[i] | mu[i], sigma));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + lognormal_lpdf(day[i] | mu[i], sigma);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ lognormal( mu[i], sigma );
}
"
m1.10 <- stan(model_code=stanmodel1.10, data=d)
SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000125 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.25 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 3.14968 seconds (Warm-up)
Chain 1: 2.42543 seconds (Sampling)
Chain 1: 5.57511 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.000141 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.41 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 2.67482 seconds (Warm-up)
Chain 2: 1.71738 seconds (Sampling)
Chain 2: 4.3922 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000107 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.07 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 2.77018 seconds (Warm-up)
Chain 3: 1.81761 seconds (Sampling)
Chain 3: 4.58779 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000139 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.39 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 3.14339 seconds (Warm-up)
Chain 4: 3.44874 seconds (Sampling)
Chain 4: 6.59213 seconds (Total)
Chain 4: There were 35 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problemsprecis(m1.10, depth = 2)shinystan::launch_shinystan(m1.10)post <- as.data.frame(m1.10)
head(post)now need to make predictions
post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu <- post_logavgdays %>%
summarise_all(mean)
mu_sigma <- mean(post$sigma)post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_pposterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
mu=t(mu),
p=t(mu_p),
sigma=mu_sigma)plot, include original data
post_plot <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=plnorm(day, meanlog=mu, sdlog=sigma) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAweek2
mv norm for rate and p?
Look at Stancode from chapter 14 mvnorm
rethinking::stancode(ulam(
alist(
wait ~ normal( mu , sigma ),
mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
c(a_cafe,b_cafe)[cafe] ~ multi_normal( c(a,b) , Rho , sigma_cafe ),
a ~ normal(5,2),
b ~ normal(-1,0.5),
sigma_cafe ~ exponential(1),
sigma ~ exponential(1),
Rho ~ lkj_corr(2)
) , data=d , sample = FALSE ))
data{
vector[200] wait;
int afternoon[200];
int cafe[200];
}
parameters{
vector[20] b_cafe;
vector[20] a_cafe;
real a;
real b;
vector<lower=0>[2] sigma_cafe;
real<lower=0> sigma;
corr_matrix[2] Rho;
}
model{
vector[200] mu;
Rho ~ lkj_corr( 2 );
sigma ~ exponential( 1 );
sigma_cafe ~ exponential( 1 );
b ~ normal( -1 , 0.5 );
a ~ normal( 5 , 2 );
{
vector[2] YY[20];
vector[2] MU;
MU = [ a , b ]';
for ( j in 1:20 ) YY[j] = [ a_cafe[j] , b_cafe[j] ]';
YY ~ multi_normal( MU , quad_form_diag(Rho , sigma_cafe) );
}
for ( i in 1:200 ) {
mu[i] = a_cafe[cafe[i]] + b_cafe[cafe[i]] * afternoon[i];
}
wait ~ normal( mu , sigma );
}Make shape exponential
d <- list(N=nrow(germ.stdi),
germ=germ.stdi$germ,
temps=as.numeric(as.factor(germ.stdi$temps)),
day=germ.stdi$day)
stanmodel1.7mv <-
"
data{
int<lower=1> N; // number of observations
int germ[N];
vector[N] day;
int temps[N];
}
parameters{
real<lower=0> shape; // should set lower bound
vector[8] a_temp; // alpha for the gamma curve, one for each temp
vector[8] ap_temp; // alpha for the proportion dormant, one for each temp
real a;
real ap;
vector<lower=0>[2] sigma_temp;
corr_matrix[2] Rho;
}
model{
vector[N] p;
vector[N] mu;
Rho ~ lkj_corr(2);
sigma_temp ~ exponential(1);
a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
ap ~ normal( 0, 1.5 );
shape ~ exponential(.5);
{
vector[2] YY[8];
vector[2] MU;
MU = [a, ap]';
for (j in 1:8) YY[j] = [a_temp[j], ap_temp[j]]';
YY ~ multi_normal(MU, quad_form_diag(Rho, sigma_temp));
}
for (i in 1:N) {
p[i] = ap_temp[temps[i]];
p[i] = inv_logit(p[i]);
}
for ( i in 1:N ) {
mu[i] = a_temp[temps[i]];
// apply the inverse link function
mu[i] = shape * exp(-(mu[i]));
}
for ( i in 1:N )
if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
for ( i in 1:N )
if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
for ( i in 1:N )
if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7mv <- stan(model_code=stanmodel1.7mv, data=d, chains=4, cores=4, control=list(adapt_delta=.99))starting worker pid=95107 on localhost:11147 at 12:45:01.786
starting worker pid=95122 on localhost:11147 at 12:45:02.133
starting worker pid=95137 on localhost:11147 at 12:45:02.451
starting worker pid=95152 on localhost:11147 at 12:45:02.774
SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000776 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.001002 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000718 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000593 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.93 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 27.0793 seconds (Warm-up)
Chain 1: 26.5605 seconds (Sampling)
Chain 1: 53.6397 seconds (Total)
Chain 1:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 27.0663 seconds (Warm-up)
Chain 3: 27.5972 seconds (Sampling)
Chain 3: 54.6635 seconds (Total)
Chain 3:
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 27.7892 seconds (Warm-up)
Chain 2: 28.2471 seconds (Sampling)
Chain 2: 56.0363 seconds (Total)
Chain 2:
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 34.5539 seconds (Warm-up)
Chain 4: 22.4221 seconds (Sampling)
Chain 4: 56.976 seconds (Total)
Chain 4: There were 100 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems100 divergent transitions…
precis(m1.7mv, depth = 3)post <- as.data.frame(m1.7mv)
post_logavgdays <- post %>% select(shape, starts_with("a_"
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>%
mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_ratemu_shape <- mean(post$shape)
post_logitp <- post %>% select(starts_with("ap_")) #logit p dormant
mu_p <- post_logitp %>%
summarize_all(mean) %>%
mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
rate=t(mu_rate),
p=t(mu_p))
post_plot1.7mv <- expand_grid(posterior_coef, day=1:28) %>%
mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))
stdi.plot <- germ %>% filter(pops=="STDI") %>%
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop_germ=cumulative_germ/total_seeds)
post_plot1.7mv %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
geom_line() +
geom_point(data=stdi.plot)
NAmake a comparison plot
post_plot_compare <- bind_rows(old=post_plot1.7, new=post_plot1.7mv, .id = "fit")
post_plot_compare %>%
ggplot(aes(x=day,y=prop_germ,color=temps,group=str_c(temps,fit),linetype=fit)) +
geom_line() +
geom_point(data=stdi.plot, aes(x=day,y=prop_germ,color=temps), inherit.aes = FALSE)
NAoverall not as good
extend to multiple species/pops
modify model 1.7a
length(unique(germone$pops))
[1] 24how can we switch this around to use day as a predictor?
LS0tCnRpdGxlOiAiZ2VybWluYXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyfQpsaWJyYXJ5KHJldGhpbmtpbmcpCmxpYnJhcnkoYnJtcykKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKYGBge3J9Cmdlcm0gPC0gcmVhZF9jc3YoIi4uL0RpbWVuc2lvbnMvaHlkcm90aGVybWFsLWFsbC1zcGVjaWVzL2RhdGEvbGlnaHRfcm91bmQxX3RhbGwuY3N2IikgJT4lCiAgZmlsdGVyKHdwcyA9PSAwKSAlPiUKICBzZWxlY3QocG9wcywgdGVtcHMsIHRvdGFsX3NlZWRzLCBnZXJtLCBkYXksIGN1bXVsYXRpdmVfZ2VybSkKZ2VybQpgYGAKCndoYXQgaWYgbmVlZCBvbmUgZXZlbnQgcGVyIHJvdzoKCmBgYHtyfQpvbmVfcGVyX3JvdyA8LSBmdW5jdGlvbihkZikgewogIHRvdGFsX3NlZWQgPC0gbWF4KGRmJHRvdGFsX3NlZWRzLCBzdW0oZGYkZ2VybSkpCiAgbmV3ZGF0YSA8LSB0aWJibGUoaWQ9MTp0b3RhbF9zZWVkLCBnZXJtPTAsIGRheT1tYXgoZGYkZGF5KSkKICBkZiA8LSBkZiAlPiUgZmlsdGVyKGdlcm0+MCkKICBjb3VudCA8LSAxCiAgaWYgKG5yb3coZGYpID4gMCkgewogICAgZm9yIChpIGluIDE6bnJvdyhkZikpIHsgIyB3ZSBsb29rIGF0IGVhY2ggcm93IG9mIHRoZSBkZiB3aGVyZSBnZXJtaW5hdGlvbiBvY2N1cmVkCiAgICAgIGZvciAoaiBpbiAxOmRmJGdlcm1baV0pIHsgIyBub3cgdXBkYXRlIHRoZSBuZXdkYXRhIHRvIHJlZmxlY3QgdGhlIGdlcm1pYW50aW9uIG9mIGVhY2ggc2VlZAogICAgICAgIG5ld2RhdGEkZ2VybVtjb3VudF0gPC0gMQogICAgICAgIG5ld2RhdGEkZGF5W2NvdW50XT1kZiRkYXlbaV0KICAgICAgICBjb3VudCA8LSBjb3VudCsxICMgY291bnQga2VlcHMgdHJhY2sgb2Ygd2hpY2ggaW5kaXZpZHVhbCB3ZSBhcmUgYXQgaW4gdGhlIG5ldyBkYXRhCiAgICAgIH0gIyBmb3IgagogICAgfSAjIGZvciBpCiAgfSAjIGlmIAogIHJldHVybihuZXdkYXRhKQp9CgpnZXJtb25lIDwtIGdlcm0gJT4lIGdyb3VwX2J5KHBvcHMsIHRlbXBzKSAlPiUKICBzZWxlY3QoLWN1bXVsYXRpdmVfZ2VybSkgJT4lICMgbm90IG5lZWRlZCBpbiB0aGlzIGVuY29kaW5nIChJIHRoaW5rLi4uaW4gYW55IGNhc2Ugd291bGQgbmVlZCB0byBiZSByZWNhbGN1bGF0ZWQpCiAgbmVzdCgpICU+JQogIG11dGF0ZShuZXdkYXRhPW1hcChkYXRhLCBvbmVfcGVyX3JvdykpICU+JQogIHNlbGVjdCgtZGF0YSkgJT4lCiAgdW5uZXN0KG5ld2RhdGEpCgpnZXJtb25lCmBgYAoKIyMgMQoKSXMgZWZmZWN0IG9mIHRlbXAgfmxpbmVhciBvbiBjdW0gZ2VybT8KCmBgYHtyfQpnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiLCBkYXk9PTI4KSAlPiUKICBnZ3Bsb3QoYWVzKHg9dGVtcHMseT1jdW11bGF0aXZlX2dlcm0pKSArCiAgZ2VvbV9jb2woKQpgYGAKbm8sIHNvIHdlIGNhbid0IHRyZWF0IHRlbXAgYXMgYSBjb250aW51b3VzIGxpbmVhciBwcmVkaWN0b3IKCmBgYHtyfQpnZXJtLnN0ZGkgPC0gZ2VybW9uZSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIGRwbHlyOjpzZWxlY3QoLXBvcHMpCmdlcm0uc3RkaQpgYGAKCiMjIDIKCiMjIyBjZW5zb3JlZCBleHBvbmVudGlhbAoKIyMjIyByZXRoaW5raW5nCgpleHBvbmVudGlhbCByYXRlIGN1cnZlLCBjZW5zb3Jpbmcgc2VlZCB0aGF0IGRvbid0IGdlcm1pbmF0ZS4KYGBge3J9CmQgPC0gbGlzdChnZXJtPWdlcm0uc3RkaSRnZXJtLCAKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgptMS4xIDwtIHVsYW0oCiAgYWxpc3QoCiAgICBkYXkgfCBnZXJtPT0xIH4gZXhwb25lbnRpYWwoIGxhbWJkYSksCiAgICBkYXkgfCBnZXJtPT0wIH4gY3VzdG9tKGV4cG9uZW50aWFsX2xjY2RmKCAhWSB8IGxhbWJkYSkpLAogICAgbGFtYmRhIDwtIDEuMCAvIG11LAogICAgbG9nKG11KSA8LSBhW3RlbXBzXSwKICAgIGFbdGVtcHNdIH4gbm9ybWFsKDAsMSkpLAogIGRhdGE9ZCwKICBjaGFpbnM9NCwKICBjb3JlcyA9IDQKKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuMSwgZGVwdGggPSAyKQpgYGAKVGhlIGFib3ZlIHJlcHJlc2VudCBsb2cobWVhbiB0aW1lIHRvIGdlcm1pbmF0aW9uKQoKYGBge3J9CmV4cCg1LjQ1KQpleHAoMi40OCkKYGBgCgoKcG9zdGVyaW9yCmBgYHtyfQpwcmVkZGF0YSA8LSBleHBhbmRfZ3JpZCh0ZW1wcz0xOjgsIGRheT0xOjI4KQpwcmVkIDwtIGxpbmsobTEuMSwgZGF0YSA9IHByZWRkYXRhKQpzdHIocHJlZCkKYGBgCgptdSBpcyBhdmVyYWdlIGRheXMgdG8gZ2VybSwgbGFtYmRhIGlzIHJhdGUgcGFyYW1ldGVyLiAgbmVpdGhlciBjaGFuZ2Ugb3ZlciB0aW1lLCBvZiBjb3Vyc2UsIHNvIGhhdmluZyBkYXlzIGRvZXNuJ3QgcmVhbGx5IG1ha2Ugc2Vuc2UuCgpgYGB7cn0KcHJlZGRhdGEkbXUgPC0gYXBwbHkocHJlZCRtdSwyLG1lYW4pCnByZWRkYXRhJGxvdyA8LSBhcHBseShwcmVkJG11LDIsSFBESSlbMSxdCnByZWRkYXRhJGhpZ2ggPC0gYXBwbHkocHJlZCRtdSwgMiwgSFBESSlbMixdCmBgYAoKU2luZ2xlIHRlbXAuICBEb24ndCBuZWVkIGRheQpwb3N0ZXJpb3IKYGBge3J9CnByZWRkYXRhIDwtIGV4cGFuZF9ncmlkKHRlbXBzPTMpCnByZWQgPC0gbGluayhtMS4xLCBkYXRhID0gcHJlZGRhdGEpCnN0cihwcmVkKQpgYGAKCmhvdyB0byBjb252ZXJ0IHRvIHByb2JzPyB1c2UgcGV4cC4KCmBgYHtyfQpwcmVkcHJvYnMgPC0gcGV4cCgxOjI4LHJhdGU9cHJlZCRsYW1iZGFbMV0pCmBgYAoKCmBgYHtyfQojcGxvdCh4PTE6MjgseT1wcmVkcHJvYnMsIHR5cGU9ImwiKSAjIGNyYXNoZXMgb24gcGxvdApgYGAKCmV2ZW4gdGhvdWdoIGl0IGlzbid0IHVzaW5nIGRheSwgaW5jbHVkaW5nIGRheSBpbiB0aGUgcHJlZGljdGlvbiBkYXRhIGZyYW1lIHdpbGwgaGVscCBtZSBrZWVwIGRhdGEgaW4gdGhlIGNvcnJlY3QgZm9ybWF0LiAgTWF5YmUuCgpgYGB7cn0KcHJlZGRhdGEgPC0gZXhwYW5kX2dyaWQodGVtcHM9MTo4LCBkYXk9MToyOCkKcHJlZCA8LSBsaW5rKG0xLjEsIGRhdGEgPSBwcmVkZGF0YSkKc3RyKHByZWQpCmBgYAoKYGBge3J9CnByZWRyZXN1bHRzIDwtIHByZWRkYXRhICU+JQogIG11dGF0ZShsYW1iZGE9YXMubGlzdChhcy5kYXRhLmZyYW1lKHByZWQkbGFtYmRhKSkpCnByZWRyZXN1bHRzCmBgYAoKYGBge3J9CnByZWRyZXN1bHRzIDwtIHByZWRyZXN1bHRzICU+JQogIG11dGF0ZShwcm9icz1tYXAyKGRheSwgbGFtYmRhLCB+IHBleHAoLngsIC55KSksCiAgICAgICAgIG11PW1hcF9kYmwocHJvYnMsIG1lYW4pLAogICAgICAgICBsb3dlcj1tYXBfZGJsKHByb2JzLCB+IEhQREkoLilbMV0gJT4lIHVubGlzdCgpKSwKICAgICAgICAgdXBwZXI9bWFwX2RibChwcm9icywgfiBIUERJKC4pWzJdKSAlPiUgdW5saXN0KCkpCnByZWRyZXN1bHRzCmBgYAoKCmBgYHtyfQpwcmVkcmVzdWx0cyAlPiUgc2VsZWN0KC1sYW1iZGEsIC1wcm9icykgJT4lCiAgbXV0YXRlKHRlbXBzPWZhY3Rvcih0ZW1wcywgbGFiZWxzPWFzLmNoYXJhY3Rlcihzb3J0KHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSkpKSAlPiUKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9bXUsY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgCmBgYAoKQWRkIHJlYWxkYXRhOgoKYGBge3J9CnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcHJlZHJlc3VsdHMgJT4lIHNlbGVjdCgtbGFtYmRhLCAtcHJvYnMpICU+JQogIG11dGF0ZSh0ZW1wcz1mYWN0b3IodGVtcHMsIGxhYmVscz1hcy5jaGFyYWN0ZXIoc29ydCh1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSkpKSkgJT4lCiAgZ2dwbG90KGFlcyh4PWRheSx5PW11LGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGFlcyh5PXByb3BfZ2VybSksIGRhdGE9c3RkaS5wbG90KQpgYGAKClBvb3IhCgojIyMgYnJtcyBjZW5zb3JpbmcgbW9kZWwKCm5lZWQgdG8gc2V0IHVwIGluZGljYXRvciBmb3IgY2Vuc29yaW5nLgpgYGB7cn0KZ2VybS5zdGRpIDwtIGdlcm0uc3RkaSAlPiUKICBtdXRhdGUoY2Vucz1pZmVsc2UoZ2VybT09MCwgInJpZ2h0IiwgIm5vbmUiKSwKICAgICAgICAgdGVtcHNjPWFzLmNoYXJhY3Rlcih0ZW1wcykgJT4lIHN0cl9wYWQod2lkdGg9MiwgcGFkPSIwIikpCmdlcm0uc3RkaQpgYGAKCmBgYHtyfQpnZXRfcHJpb3IoZGF5IHwgY2VucyhjZW5zKSB+IDAgKyB0ZW1wc2MsIGZhbWlseSA9IGV4cG9uZW50aWFsLCBkYXRhPWdlcm0uc3RkaSkKYGBgCgoKYGBge3J9Cm0xLjIgPC0gYnJtKGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgICBmYW1pbHkgPSBleHBvbmVudGlhbCgpLAogICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiwgY2xhc3M9ImIiKSwKICAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSwgKQpgYGAKCmBgYHtyfQpzdW1tYXJ5KG0xLjIpCmBgYAoKCmBgYHtyfQpwcmVkaWN0KG0xLjIpICU+JSBoZWFkCmBgYApsb2cobWVhbiB0aW1lIHRvIGdlcm1pbmF0aW9uKQoKYGBge3J9CnByZWRpY3QobTEuMiwgbmV3ZGF0YSA9IGV4cGFuZF9ncmlkKHRlbXBzYz11bmlxdWUoZ2VybS5zdGRpJHRlbXBzYyksIGNlbnM9Yygibm9uZSIsICJyaWdodCIpKSkKYGBgCgpgYGB7cn0KcHJlY2lzKG0xLjEsIGRlcHRoID0gMikKYGBgCmFnYWluLCB0aGVzZSBhcmUgbG9nKG1lYW4gdGltZSB0byBnZXJtaW5hdGlvbikKCmxhbWJkYSBzaG91bGQgYmUgMS9tdSwgc28KCmBgYHtyfQojcGxvdCgxOjI4LHBleHAoMToyOCwgMS9leHAoMi40OCkpLCB0eXBlPSJsIiwgY29sPSJyZWQiKQojbGluZXMoMToyOCxwZXhwKDE6MjgsIDEvZXhwKDUuNjApKSwgdHlwZT0ibCIsIGNvbD0iYmx1ZSIpCiNjcmFzaGVzCmBgYAoKCmBgYHtyfQpnZXRfcHJpb3IoIGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgIGZhbWlseSA9IEdhbW1hKCksCiAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSkKYGBgCgoKYGBge3J9Cm0xLjMgPC0gYnJtKGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgICBmYW1pbHkgPSBHYW1tYSgpLAogICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiwgY2xhc3M9ImIiLCBsYj0wKSwKICAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShtMS4zKQpgYGAKCmBgYHtyfQpuZXdkYXRhIDwtIGV4cGFuZF9ncmlkKHRlbXBzYz11bmlxdWUoZ2VybS5zdGRpJHRlbXBzYyksIGNlbnM9Yygibm9uZSIsICJyaWdodCIpKQpjYmluZChuZXdkYXRhLCBwcmVkaWN0KG0xLjMsIG5ld2RhdGEgPSBuZXdkYXRhKSkKYGBgCgpUaGUgcHJvYmxlbSB3aXRoIHRoZXNlIG1vZGVscyBpcyB0aGF0IHRoZXkgYXNzdW1lIHRoYXQgdGhlIGNlbnNvcmVkIG9ic2VydmF0aW9ucyBhcmUgY29udHJvbGxlZCBieSB0aGUgc2FtZSBwcm9jZXNzIHRoYXQgaXMgY29udHJvbGxpbmcgdGhlIG9uZXMgdGhhdCBkbyBnZXJtaW5hdGUuICBJIGRvIG5vdCB0aGluayB0aGF0IGlzIGNvcnJlY3QsIHdlIG5lZWQgYSBtaXh0dXJlIG1vZGVsLgoKIyMgVHJ5aW5nIGEgbWl4dHVyZSBtb2RlbCBmb3IgYSBaSSBleHBvbmVudGlhbAoKbG9vayBhdCB1bGFtIGNlbnNvcmVkIGV4cG9uZW50aWFsCmBgYHtyfQpyZXRoaW5raW5nOjpzdGFuY29kZShtMS4xKQpgYGAKCmxvb2sgYXQgYnJtcyBjZW5zb3JlZCBleHBvbmVudGlhbApgYGB7cn0Kc3RhbmNvZGUobTEuMikKYGBgCgpaZXJvIGluZmxhdGVkLCBhcyBnZW5lcmF0ZWQgYnkgdWxhbToKYGBge3J9CiAKIyBkZWZpbmUgcGFyYW1ldGVycwpwcm9iX2RyaW5rIDwtIDAuMiAjIDIwJSBvZiBkYXlzCnJhdGVfd29yayA8LSAxICMgYXZlcmFnZSAxIG1hbnVzY3JpcHQgcGVyIGRheQojIHNhbXBsZSBvbmUgeWVhciBvZiBwcm9kdWN0aW9uIApOIDwtIDM2NQojIHNpbXVsYXRlIGRheXMgbW9ua3MgZHJpbmsgCnNldC5zZWVkKDM2NSkKZHJpbmsgPC0gcmJpbm9tKCBOICwgMSAsIHByb2JfZHJpbmsgKQojIHNpbXVsYXRlIG1hbnVzY3JpcHRzIGNvbXBsZXRlZAp5IDwtICgxLWRyaW5rKSpycG9pcyggTiAsIHJhdGVfd29yayApCgpyZXRoaW5raW5nOjpzdGFuY29kZSh1bGFtKCBhbGlzdCgKICB5fHk+MCB+IGN1c3RvbSggbG9nMW0ocCkgKyBwb2lzc29uX2xwbWYoeXxsYW1iZGEpICksCiAgeXx5PT0wIH4gY3VzdG9tKCBsb2dfbWl4KCBwICwgMCAsIHBvaXNzb25fbHBtZigwfGxhbWJkYSkgKSApLAogIGxvZ2l0KHApIDwtIGFwLAogIGxvZyhsYW1iZGEpIDwtIGFsLAogIGFwIH4gZG5vcm0oLTEuNSwxKSwKICBhbCB+IGRub3JtKDEsMC41KSksCiAgZGF0YT1saXN0KHk9YXMuaW50ZWdlcih5KSkgLCBzYW1wbGU9RkFMU0UpKQpgYGAKClplcm8gaW5mbGF0ZWQgZXhwb25lbnRpYWwgbW9kZWwsIGF0dGVtcHQgMQpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjQgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhcyBmb3IgdGhlIGxpbmVhciBtb2RlbAogICAgcmVhbCBhcDsgLy8gcHJvYmFiaWxpdHkgb2YgZHJpbmsKfQptb2RlbHsKICAgIHJlYWwgcDsKICAgIHZlY3RvcltOXSBsYW1iZGE7CiAgICB2ZWN0b3JbTl0gbXU7CiAgICBhIH4gbm9ybWFsKCAwICwgMSApOwogICAgYXAgfiBub3JtYWwoIDAsIDEgKTsKICAgIHAgPSBhcDsKICAgIHAgPSBpbnZfbG9naXQocCk7IC8vIGludmVyc2UgbGluayBmdW5jdGlvbiBmb3IgcHJvYmFiaWxpdHkgZHJpbmsKICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07IC8vIGxpbmVhciBtb2RlbCBmb3IgbXUKICAgICAgICBtdVtpXSA9IGV4cChtdVtpXSk7ICAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24gcGFydCAxCiAgICAgICAgbGFtYmRhW2ldID0gMS9tdVtpXTsgLy8gaW52ZXJzZSBsaW5rIGZ1bmN0aW9uIHBhcnQgMgogICAgfQoKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwLCAwLCBleHBvbmVudGlhbF9sY2NkZihkYXlbaV0gfCBsYW1iZGFbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocCkgKyAgZXhwb25lbnRpYWxfbHBkZihkYXlbaV0gfCBsYW1iZGFbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGV4cG9uZW50aWFsKCBsYW1iZGFbaV0gKTsKfQoiCm0xLjQgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuNCwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuNCwgZGVwdGggPSAyKQpgYGAKTW9kZWwgMS40IGZpdHMgYSBzaW5nbGUgZG9ybWFuY3kgcmF0ZSBmb3IgYWxsLCBidXQgcHJvYmFibHkgd2UgbmVlZCB0byBoYXZlIHNlcGFyYXRlIGRvcm1hbmN5IHJhdGUgZm9yIGVhY2ggdGVtcAoKdGFrZSBhIGxvb2sgYXQgc3RhbmNvZGUgZm9yIGEgYmlub21pYWwgbW9kZWwgdG8gc2VlIHdoYXQgdGhpcyBsb29rcyBsaWtlCgpgYGB7cn0KZGIgPC0gbGlzdChjdW11bGF0aXZlX2dlcm09Z2VybSRjdW11bGF0aXZlX2dlcm0sIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0kdGVtcHMpKSkKCnJldGhpbmtpbmc6OnN0YW5jb2RlKHVsYW0oYWxpc3QoCiAgY3VtdWxhdGl2ZV9nZXJtIH4gZGJlcm4obXUpLAogIG11IDwtIGxvZ2l0KGEgKyBiW3RlbXBzXSksCiAgYSB+IGRub3JtKDAsMSksCiAgYlt0ZW1wc10gfiBkbm9ybSgwLDEpKSwKICBkYXRhID0gZGIsIHNhbXBsZT1GQUxTRSwgbG9nX2xpayA9IFRSVUUpKQpgYGAKVGhlcmUgaXMgbm8gaW52ZXJzZSBsb2dpdCBiZWNhdXNlIGl0IGlzIGJ1aWx0IGluIHRvIHRoZSBiZXJub3VsbGkgZnVuY3Rpb24KClRyeSB0byB3cml0ZSBhIG1vZGVsIHdpdGggZGlmZmVyZW50IGRvcm1hbmN5IHByb2JzIGZvciB0aGUgZGlmZmVyZW50IHRlbXBzLgpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjUgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZXhwb25lbnRpYWwgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbGFtYmRhOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIDEgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAxICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgbXVbaV0gPSBleHAobXVbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIGxhbWJkYVtpXSA9IDEvbXVbaV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZXhwb25lbnRpYWxfbGNjZGYoZGF5W2ldIHwgbGFtYmRhW2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgIGV4cG9uZW50aWFsX2xwZGYoZGF5W2ldIHwgbGFtYmRhW2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBleHBvbmVudGlhbCggbGFtYmRhW2ldICk7Cn0KIgptMS41IDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjUsIGRhdGE9ZCkKYGBgCgpgYGB7cn0KcHJlY2lzKG0xLjUsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjUpCmBgYAoKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuNSkKaGVhZChwb3N0KQpgYGAKCmBgYHtyfQpwb3N0LmNvciA8LSBjb3IocG9zdCkgJT4lIHJvdW5kKDIpCnBvc3QuY29yW2Ficyhwb3N0LmNvcik8LjI1XSA8LSBOQQpwb3N0LmNvcgpgYGAKCgpub3cgbmVlZCB0byBtYWtlIHByZWRpY3Rpb25zCgpgYGB7cn0KcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9sYW1iZGEgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pICU+JSAjIHBvc3RlcmlvciBtZWFuCiAgbXV0YXRlX2FsbCh+IDEvZXhwKC4pKSAjIGNvbnZlcnQgdG8gbGFtYmRhIChpbnZlcnNlIGxpbmspCm11X2xhbWJkYQpgYGAKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZSh0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICBsYW1iZGE9dChtdV9sYW1iZGEpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCmBgYAoKcGxvdCwgaW5jbHVkZSBvcmlnaW5hbCBkYXRhCmBgYHtyfQpwb3N0X3Bsb3QgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBleHAoZGF5LCByYXRlPWxhbWJkYSkgKiAoMS1wKSkKCnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcG9zdF9wbG90ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCmFkanVzdCBwcmlvcnMKCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuNWEgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZXhwb25lbnRpYWwgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbGFtYmRhOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIDIgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAyICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgbXVbaV0gPSBleHAobXVbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIGxhbWJkYVtpXSA9IDEvbXVbaV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZXhwb25lbnRpYWxfbGNjZGYoZGF5W2ldIHwgbGFtYmRhW2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgIGV4cG9uZW50aWFsX2xwZGYoZGF5W2ldIHwgbGFtYmRhW2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBleHBvbmVudGlhbCggbGFtYmRhW2ldICk7Cn0KIgptMS41YSA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS41YSwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuNWEsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjVhKQpgYGAKCgpgYGB7cn0KcG9zdCA8LSBhcy5kYXRhLmZyYW1lKG0xLjVhKQoKcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9sYW1iZGEgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pICU+JSAjIHBvc3RlcmlvciBtZWFuCiAgbXV0YXRlX2FsbCh+IDEvZXhwKC4pKSAjIGNvbnZlcnQgdG8gbGFtYmRhIChpbnZlcnNlIGxpbmspCm11X2xhbWJkYQoKcG9zdF9sb2dpdHAgPC0gcG9zdCAlPiUgc2VsZWN0KHN0YXJ0c193aXRoKCJhcCIpKSAjbG9naXQgcCBkb3JtYW50Cm11X3AgPC0gcG9zdF9sb2dpdHAgJT4lCiAgc3VtbWFyaXplX2FsbChtZWFuKSAlPiUKICBtdXRhdGVfYWxsKGludl9sb2dpdCkKbXVfcAoKcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIGxhbWJkYT10KG11X2xhbWJkYSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKCnBvc3RfcGxvdCA8LSBleHBhbmRfZ3JpZChwb3N0ZXJpb3JfY29lZiwgZGF5PTE6MjgpICU+JQogIG11dGF0ZShwcm9wX2dlcm09cGV4cChkYXksIHJhdGU9bGFtYmRhKSAqICgxLXApKQoKc3RkaS5wbG90IDwtIGdlcm0gJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSAKICBzZWxlY3QoZGF5LCB0ZW1wcywgY3VtdWxhdGl2ZV9nZXJtLCB0b3RhbF9zZWVkcykgJT4lCiAgbXV0YXRlKHRlbXBzPWFzLmZhY3Rvcih0ZW1wcyksCiAgICAgICAgIHByb3BfZ2VybT1jdW11bGF0aXZlX2dlcm0vdG90YWxfc2VlZHMpCgpwb3N0X3Bsb3QgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYApzdGlsbCBub3QgZ3JlYXQKCiMjIyBnYW1tYQoKTG9vayBhdCBnYW1tYSBzdGFuY29kZSBmcm9tIGJybXM6CgpgYGB7cn0KbWFrZV9zdGFuY29kZShkYXkgfiB0ZW1wc2MsIAogICAgICAgICAgICAgZmFtaWx5PSJHYW1tYSIsCiAgICAgICAgICAgICBwcmlvcj1zZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiksCiAgICAgICAgICAgICBkYXRhPWdlcm0uc3RkaSwKICAgICAgICAgICAgIHNhbXBsZT1GQUxTRSkKYGBgCgoKVHJ5IGRldmVsb3BpbmcgYSBaSSBnYW1tYSBtb2RlbDoKCnN0YW4gdXNlcyB0aGF0IGdhbW1hKGFscGhhLCBiZXRhKSBwYXJhbWV0ZXJpemF0aW9uIHdoZXJlICJhbHBoYSIgaXMgc2hhcGUiIGFuZCAiYmV0YSIgaXMgInJhdGUiLiAgU28gdHJ5IGhhdmluZyB0aGUgc2hhcGUgYmUgdGhlIHNhbWUgZm9yIGFsbCB0ZW1wcy4KYGBge3J9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS43IDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC41ICk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICBzaGFwZSB+IG5vcm1hbCgwLCAuNSk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBtdVtpXSA9IGFbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9uCiAgICAgICAgbXVbaV0gPSBzaGFwZSAqIGV4cCgtKG11W2ldKSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwW2ldLCAwLCBnYW1tYV9sY2NkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGdhbW1hKCBzaGFwZSwgbXVbaV0gKTsKfQoiCm0xLjcgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuNywgZGF0YT1kLCBjaGFpbnM9NCwgY29yZXM9NCwgY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPS45OSkpCmBgYCAKCgpgYGB7cn0KcHJlY2lzKG0xLjcsIGRlcHRoID0gMikKYGBgCgoKYGBge3J9CnNoaW55c3Rhbjo6bGF1bmNoX3NoaW55c3RhbihtMS43KQpgYGAKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuNykKaGVhZChwb3N0KQpgYGAKCgpub3cgbmVlZCB0byBtYWtlIHByZWRpY3Rpb25zCgpgYGB7cn0KcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzaGFwZSwgc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCm11X3NoYXBlCmBgYAoKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKCnBsb3QgcHJlZGljdGlvbnMgYW5kIG9yaWdpbmFsIGRhdGEKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZShzaGFwZT1tdV9zaGFwZSwgdGVtcHM9YXMuZmFjdG9yKHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgICAgICAgICAgICAgICAgcmF0ZT10KG11X3JhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCgpwb3N0X3Bsb3QxLjcgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuNyAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCgojIyMjIHBsYXlpbmcgd2l0aCBwcmlvcnM6CgpNYWtlIHNoYXBlIGV4cG9uZW50aWFsCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuN2EgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgcmVhbDxsb3dlcj0wPiBzaGFwZTsgLy8gc2hvdWxkIHNldCBsb3dlciBib3VuZAogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZ2FtbWEgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbXU7CiAgICBhIH4gbm9ybWFsKCAwICwgLjUgKTsgLy9uYXJyb3cgcHJpb3JzIHRvIG92ZXJjb21lIGRpdmVyZ2VudCB0cmFuc2l0aW9ucwogICAgYXAgfiBub3JtYWwoIDAsIDEuNSApOyAKICAgIHNoYXBlIH4gZXhwb25lbnRpYWwoLjUpOyAKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBtdVtpXSA9IGFbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9uCiAgICAgICAgbXVbaV0gPSBzaGFwZSAqIGV4cCgtKG11W2ldKSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwW2ldLCAwLCBnYW1tYV9sY2NkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGdhbW1hKCBzaGFwZSwgbXVbaV0gKTsKfQoiCm0xLjdhIDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjdhLCBkYXRhPWQsIGNoYWlucz00LCBjb3Jlcz00LCBjb250cm9sPWxpc3QoYWRhcHRfZGVsdGE9Ljk5KSkKYGBgIAoKCmBgYHtyfQpwcmVjaXMobTEuN2EsIGRlcHRoID0gMikKYGBgCgoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS43YSkKCnBvc3RfbG9nYXZnZGF5cyA8LSBwb3N0ICU+JSBzZWxlY3Qoc2hhcGUsIHN0YXJ0c193aXRoKCJhWyIKKSkgIyB0aGVzZSBhcmUgdGhlIHNoYXBlIGFuZCB0aGUgbG9nKG1lYW4oYXZnIHRpbWUgdG8gZ2VybWluYXRlKSkKbXVfcmF0ZSA8LSBwb3N0X2xvZ2F2Z2RheXMgJT4lIAogIG11dGF0ZShhY3Jvc3MoLXNoYXBlLCB+IHNoYXBlKmV4cCgtKC4pKSkgKSAlPiUKICBzdW1tYXJpemUoYWNyb3NzKGV2ZXJ5dGhpbmcoKSwgbWVhbikpICU+JSBzZWxlY3QoLXNoYXBlKSAjIHBvc3RlcmlvciBtZWFuCm11X3JhdGUKbXVfc2hhcGUgPC0gbWVhbihwb3N0JHNoYXBlKQoKcG9zdF9sb2dpdHAgPC0gcG9zdCAlPiUgc2VsZWN0KHN0YXJ0c193aXRoKCJhcCIpKSAjbG9naXQgcCBkb3JtYW50Cm11X3AgPC0gcG9zdF9sb2dpdHAgJT4lCiAgc3VtbWFyaXplX2FsbChtZWFuKSAlPiUKICBtdXRhdGVfYWxsKGludl9sb2dpdCkKbXVfcAoKcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHNoYXBlPW11X3NoYXBlLCB0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICByYXRlPXQobXVfcmF0ZSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKCnBvc3RfcGxvdDEuN2EgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuN2EgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYApTaGFwZSBwYXJhbWV0ZXIgaXMgbGFyZ2VyIHRoYW4gaW4gMS43LiAgVGhlIGZpdCBpcyBzbGlnaHRseSBiZXR0ZXIgYXQgbGVhc3QgZm9yIDIwIGFuZCAyNS4KCm1ha2UgYSBjb21wYXJpc29uIHBsb3QKCmBgYHtyfQpwb3N0X3Bsb3RfY29tcGFyZSA8LSBiaW5kX3Jvd3Mob2xkPXBvc3RfcGxvdDEuNywgbmV3PXBvc3RfcGxvdDEuN2EsIC5pZCA9ICJmaXQiKSAKcG9zdF9wbG90X2NvbXBhcmUgJT4lCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD1zdHJfYyh0ZW1wcyxmaXQpLGxpbmV0eXBlPWZpdCkpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCwgYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzKSwgaW5oZXJpdC5hZXMgPSBGQUxTRSkKICAKYGBgCkp1c3QgYSBoYWlyIGJldHRlci4uLgoKY29udGludWUgdG8gd29yayB3aXRoIHByaW9ycwpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjdiIDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC43NSApOyAvLyBkaXZlcmdlbnQgd2l0aCAxIGFuZCB3aXRoIC43NQogICAgYXAgfiBub3JtYWwoIDAsIDEuNSApOyAKICAgIHNoYXBlIH4gZXhwb25lbnRpYWwoMS41KTsgLy8gZGl2ZXJnZW50IGlmIHRoaXMgaXMgLjUgYW5kIGEgaXMgMCwgMQogICAgZm9yIChpIGluIDE6TikgewogICAgICAgIHBbaV0gPSBhcFt0ZW1wc1tpXV07CiAgICAgICAgcFtpXSA9IGludl9sb2dpdChwW2ldKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgIC8vIGFwcGx5IHRoZSBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgICAgICBtdVtpXSA9IHNoYXBlICogZXhwKC0obXVbaV0pKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMCApIHRhcmdldCArPSBsb2dfbWl4KHBbaV0sIDAsIGdhbW1hX2xjY2RmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSkpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgZ2FtbWFfbHBkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgZGF5W2ldIH4gZ2FtbWEoIHNoYXBlLCBtdVtpXSApOwp9CiIKbTEuN2IgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuN2IsIGRhdGE9ZCwgY2hhaW5zPTQsIGNvcmVzPTQsIGNvbnRyb2w9bGlzdChhZGFwdF9kZWx0YT0uOTkpKQpgYGAgCgoKYGBge3J9CnByZWNpcyhtMS43YiwgZGVwdGggPSAyKQpgYGAKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuN2IpCgpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBzaGFwZSBhbmQgdGhlIGxvZyhtZWFuKGF2ZyB0aW1lIHRvIGdlcm1pbmF0ZSkpCm11X3JhdGUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBtdXRhdGUoYWNyb3NzKC1zaGFwZSwgfiBzaGFwZSpleHAoLSguKSkpICkgJT4lCiAgc3VtbWFyaXplKGFjcm9zcyhldmVyeXRoaW5nKCksIG1lYW4pKSAlPiUgc2VsZWN0KC1zaGFwZSkgIyBwb3N0ZXJpb3IgbWVhbgptdV9yYXRlCm11X3NoYXBlIDwtIG1lYW4ocG9zdCRzaGFwZSkKCnBvc3RfbG9naXRwIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYXAiKSkgI2xvZ2l0IHAgZG9ybWFudAptdV9wIDwtIHBvc3RfbG9naXRwICU+JQogIHN1bW1hcml6ZV9hbGwobWVhbikgJT4lCiAgbXV0YXRlX2FsbChpbnZfbG9naXQpCm11X3AKCnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZShzaGFwZT1tdV9zaGFwZSwgdGVtcHM9YXMuZmFjdG9yKHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgICAgICAgICAgICAgICAgcmF0ZT10KG11X3JhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCgpwb3N0X3Bsb3QgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdCAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCm5vIGltcHJvdmVtZW50IG92ZXIgMS43YQoKIyMjIGJybXMgaHVyZGxlIGdhbW1hCgpgYGB7cn0KZ2V0X3ByaW9yKGJmKGRheSB+IDAgKyB0ZW1wc2MsIGh1IH4gMCArIHRlbXBzYyksCiAgICAgICAgICBmYW1pbHkgPSBodXJkbGVfZ2FtbWEsCiAgICAgICAgICBkYXRhPWdlcm0uc3RkaSkKYGBgCgpgYGB7cn0KbWFrZV9zdGFuY29kZShiZihkYXkgfiAwICsgdGVtcHNjLCBodSB+IDAgKyB0ZW1wc2MpLAogICAgICAgICAgZmFtaWx5ID0gaHVyZGxlX2dhbW1hLAogICAgICAgICAgcHJpb3I9c2V0X3ByaW9yKCJub3JtYWwoMCwxKSIpLAogICAgICAgICAgZGF0YT1nZXJtLnN0ZGkpCmBgYAoKCmBgYHtyfQptMS44IDwtIGJybShiZihkYXkgfiAwICsgdGVtcHNjLCBodSB+IDAgKyB0ZW1wc2MpLAogICAgICAgICAgZmFtaWx5ID0gaHVyZGxlX2dhbW1hLAogICAgICAgICAgcHJpb3I9YyhzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiksCiAgICAgICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDAsMSkiLCBkcGFyPSJodSIpKSwgIyBtdXN0IHNwZWNpZnkgdGhpcyEKICAgICAgICAgIGRhdGE9Z2VybS5zdGRpLCBjb3JlcyA9IDQpCmBgYAoKCmBgYHtyfQptMS44CmBgYAoKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuOCkKaGVhZChwb3N0KQpgYGAKCndlaXJkIGZpdCEKCm5vdyBuZWVkIHRvIG1ha2UgcHJlZGljdGlvbnMKCmBgYHtyfQpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYl90ZW1wIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCmBgYAoKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImJfaHUiKSkgI2xvZ2l0IHAgZG9ybWFudAptdV9wIDwtIHBvc3RfbG9naXRwICU+JQogIHN1bW1hcml6ZV9hbGwobWVhbikgJT4lCiAgbXV0YXRlX2FsbChpbnZfbG9naXQpCm11X3AKYGBgCgpgYGB7cn0KcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHNoYXBlPW11X3NoYXBlLCB0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICByYXRlPXQobXVfcmF0ZSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKYGBgCgpwbG90IHByZWRpY3Rpb25zIGFuZCBvcmlnaW5hbCBkYXRhCmBgYHtyfQoKCnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcG9zdF9wbG90ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCiMjIGFsbG93IGJvdGggc2hhcGUgYW5kIHJhdGUgb2YgWkkgZ2FtbWEgdG8gYmUgcHJlZGljdGVkIGJ5IHRlbXAKCkRvbid0IHJ1bi4gIERvZXNuJ3Qgc2FtcGxlIHdlbGwgYmVjYXVzZSBzaGFwZSBhbmQgcmF0ZSBhcmUgY29ycmVsYXRlZC4KYGBge3IsIGV2YWw9RkFMU0V9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS45IDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHZlY3Rvcls4XSBhczsgLy8gYWxwaGEgZm9yIGdhbWEgc2hhcGUKICAgIHZlY3Rvcls4XSBhcjsgLy8gYWxwaGEgZm9yIHRoZSBnYW1tYSByYXRlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIHNoYXBlOwogICAgdmVjdG9yW05dIHJhdGU7CiAgICBhcyB+IG5vcm1hbCggMCwgMSApOwogICAgYXIgfiBub3JtYWwoIDAsIDEgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICAKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBzaGFwZVtpXSA9IGFzW3RlbXBzW2ldXTsKICAgICAgICByYXRlW2ldID0gYXJbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9ucwogICAgICAgIHNoYXBlW2ldID0gZXhwKHNoYXBlW2ldKTsKICAgICAgICByYXRlW2ldID0gZXhwKHJhdGVbaV0pOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGVbaV0sIHJhdGVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGVbaV0sIHJhdGVbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgZGF5W2ldIH4gZ2FtbWEoIHNoYXBlW2ldLCByYXRlIFtpXSApOwp9CmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKCn0KIgptMS45IDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjksIGRhdGE9ZCwgY2hhaW5zPTQsIGNvcmVzPTQsIGNvbnRyb2w9bGlzdChhZGFwdF9kZWx0YT0uOTkpKQpgYGAgCgoKYGBge3J9CnByZWNpcyhtMS45LCBkZXB0aCA9IDIpCmBgYAoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS45KQoKcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzaGFwZSwgc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCgpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCgpwb3N0ZXJpb3JfY29lZiA8LSB0aWJibGUoc2hhcGU9bXVfc2hhcGUsIHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIHJhdGU9dChtdV9yYXRlKSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApKQoKcG9zdF9wbG90IDwtIGV4cGFuZF9ncmlkKHBvc3Rlcmlvcl9jb2VmLCBkYXk9MToyOCkgJT4lCiAgbXV0YXRlKHByb3BfZ2VybT1wZ2FtbWEoZGF5LCBzaGFwZT1tdV9zaGFwZSwgcmF0ZT1yYXRlKSAqICgxLXApKQoKc3RkaS5wbG90IDwtIGdlcm0gJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSAKICBzZWxlY3QoZGF5LCB0ZW1wcywgY3VtdWxhdGl2ZV9nZXJtLCB0b3RhbF9zZWVkcykgJT4lCiAgbXV0YXRlKHRlbXBzPWFzLmZhY3Rvcih0ZW1wcyksCiAgICAgICAgIHByb3BfZ2VybT1jdW11bGF0aXZlX2dlcm0vdG90YWxfc2VlZHMpCgpwb3N0X3Bsb3QgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYAoKIyMgbG9nIG5vcm1hbAoKYGBge3J9Cm1ha2Vfc3RhbmNvZGUoZGF5IH4gdGVtcHNjLCBmYW1pbHkgPSBsb2dub3JtYWwoKSwgZGF0YT1nZXJtLnN0ZGkpCmBgYAoKCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuMTAgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgbG9nbm9ybWFsIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKICAgIHJlYWw8bG93ZXI9MD4gc2lnbWE7Cn0KbW9kZWx7CiAgICB2ZWN0b3JbTl0gcDsKICAgIHZlY3RvcltOXSBtdTsKICAgIGEgfiBub3JtYWwoIDAgLCAyICk7CiAgICBzaWdtYSB+IGV4cG9uZW50aWFsICgxKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAyICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgbG9nbm9ybWFsX2xjY2RmKGRheVtpXSB8IG11W2ldLCBzaWdtYSkpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArICBsb2dub3JtYWxfbHBkZihkYXlbaV0gfCBtdVtpXSwgc2lnbWEpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGxvZ25vcm1hbCggbXVbaV0sIHNpZ21hICk7Cn0KIgptMS4xMCA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS4xMCwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuMTAsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjEwKQpgYGAKCgpgYGB7cn0KcG9zdCA8LSBhcy5kYXRhLmZyYW1lKG0xLjEwKQpoZWFkKHBvc3QpCmBgYAoKbm93IG5lZWQgdG8gbWFrZSBwcmVkaWN0aW9ucwoKYGBge3J9CnBvc3RfbG9nYXZnZGF5cyA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgbG9nKG1lYW4oYXZnIHRpbWUgdG8gZ2VybWluYXRlKSkKbXUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pIAoKbXVfc2lnbWEgPC0gbWVhbihwb3N0JHNpZ21hKQpgYGAKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZSh0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICBtdT10KG11KSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApLAogICAgICAgICAgICAgICAgICAgICAgICBzaWdtYT1tdV9zaWdtYSkKYGBgCgpwbG90LCBpbmNsdWRlIG9yaWdpbmFsIGRhdGEKYGBge3J9CnBvc3RfcGxvdCA8LSBleHBhbmRfZ3JpZChwb3N0ZXJpb3JfY29lZiwgZGF5PTE6MjgpICU+JQogIG11dGF0ZShwcm9wX2dlcm09cGxub3JtKGRheSwgbWVhbmxvZz1tdSwgc2Rsb2c9c2lnbWEpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdCAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCgojIHdlZWsyCgojIyBtdiBub3JtIGZvciByYXRlIGFuZCBwPwoKTG9vayBhdCBTdGFuY29kZSBmcm9tIGNoYXB0ZXIgMTQgbXZub3JtCmBgYHtyfQojIyBSIGNvZGUgMTQuMQphIDwtIDMuNSAgICAgICAgICAgICMgYXZlcmFnZSBtb3JuaW5nIHdhaXQgdGltZQpiIDwtICgtMSkgICAgICAgICAgICMgYXZlcmFnZSBkaWZmZXJlbmNlIGFmdGVybm9vbiB3YWl0IHRpbWUKc2lnbWFfYSA8LSAxICAgICAgICAjIHN0ZCBkZXYgaW4gaW50ZXJjZXB0cwpzaWdtYV9iIDwtIDAuNSAgICAgICMgc3RkIGRldiBpbiBzbG9wZXMKcmhvIDwtICgtMC43KSAgICAgICAjIGNvcnJlbGF0aW9uIGJldHdlZW4gaW50ZXJjZXB0cyBhbmQgc2xvcGVzCgojIyBSIGNvZGUgMTQuMgpNdSA8LSBjKCBhICwgYiApCgojIyBSIGNvZGUgMTQuMwpjb3ZfYWIgPC0gc2lnbWFfYSpzaWdtYV9iKnJobwpTaWdtYSA8LSBtYXRyaXgoIGMoc2lnbWFfYV4yLGNvdl9hYixjb3ZfYWIsc2lnbWFfYl4yKSAsIG5jb2w9MiApCgoKIyMgUiBjb2RlIDE0LjUKc2lnbWFzIDwtIGMoc2lnbWFfYSxzaWdtYV9iKSAjIHN0YW5kYXJkIGRldmlhdGlvbnMKUmhvIDwtIG1hdHJpeCggYygxLHJobyxyaG8sMSkgLCBucm93PTIgKSAjIGNvcnJlbGF0aW9uIG1hdHJpeAoKIyBub3cgbWF0cml4IG11bHRpcGx5IHRvIGdldCBjb3ZhcmlhbmNlIG1hdHJpeApTaWdtYSA8LSBkaWFnKHNpZ21hcykgJSolIFJobyAlKiUgZGlhZyhzaWdtYXMpCgojIyBSIGNvZGUgMTQuNgpOX2NhZmVzIDwtIDIwCgojIyBSIGNvZGUgMTQuNwpzZXQuc2VlZCg1KSAjIHVzZWQgdG8gcmVwbGljYXRlIGV4YW1wbGUKdmFyeV9lZmZlY3RzIDwtIE1BU1M6Om12cm5vcm0oIE5fY2FmZXMgLCBNdSAsIFNpZ21hICkKCiMjIFIgY29kZSAxNC44CmFfY2FmZSA8LSB2YXJ5X2VmZmVjdHNbLDFdCmJfY2FmZSA8LSB2YXJ5X2VmZmVjdHNbLDJdCgojIyBSIGNvZGUgMTQuMTAKc2V0LnNlZWQoMjIpCk5fdmlzaXRzIDwtIDEwCmFmdGVybm9vbiA8LSByZXAoMDoxLE5fdmlzaXRzKk5fY2FmZXMvMikKY2FmZV9pZCA8LSByZXAoIDE6Tl9jYWZlcyAsIGVhY2g9Tl92aXNpdHMgKQptdSA8LSBhX2NhZmVbY2FmZV9pZF0gKyBiX2NhZmVbY2FmZV9pZF0qYWZ0ZXJub29uCnNpZ21hIDwtIDAuNSAgIyBzdGQgZGV2IHdpdGhpbiBjYWZlcwp3YWl0IDwtIHJub3JtKCBOX3Zpc2l0cypOX2NhZmVzICwgbXUgLCBzaWdtYSApCmQgPC0gZGF0YS5mcmFtZSggY2FmZT1jYWZlX2lkICwgYWZ0ZXJub29uPWFmdGVybm9vbiAsIHdhaXQ9d2FpdCApCgoKcmV0aGlua2luZzo6c3RhbmNvZGUodWxhbSgKICAgIGFsaXN0KAogICAgICAgIHdhaXQgfiBub3JtYWwoIG11ICwgc2lnbWEgKSwKICAgICAgICBtdSA8LSBhX2NhZmVbY2FmZV0gKyBiX2NhZmVbY2FmZV0qYWZ0ZXJub29uLAogICAgICAgIGMoYV9jYWZlLGJfY2FmZSlbY2FmZV0gfiBtdWx0aV9ub3JtYWwoIGMoYSxiKSAsIFJobyAsIHNpZ21hX2NhZmUgKSwKICAgICAgICBhIH4gbm9ybWFsKDUsMiksCiAgICAgICAgYiB+IG5vcm1hbCgtMSwwLjUpLAogICAgICAgIHNpZ21hX2NhZmUgfiBleHBvbmVudGlhbCgxKSwKICAgICAgICBzaWdtYSB+IGV4cG9uZW50aWFsKDEpLAogICAgICAgIFJobyB+IGxral9jb3JyKDIpCiAgICApICwgZGF0YT1kICwgc2FtcGxlID0gRkFMU0UgKSkKYGBgCgoKCk1ha2Ugc2hhcGUgZXhwb25lbnRpYWwKYGBge3J9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS43bXYgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgcmVhbDxsb3dlcj0wPiBzaGFwZTsgLy8gc2hvdWxkIHNldCBsb3dlciBib3VuZAogICAgdmVjdG9yWzhdIGFfdGVtcDsgLy8gYWxwaGEgZm9yIHRoZSBnYW1tYSBjdXJ2ZSwgb25lIGZvciBlYWNoIHRlbXAKICAgIHZlY3Rvcls4XSBhcF90ZW1wOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKICAgIHJlYWwgYTsKICAgIHJlYWwgYXA7CiAgICB2ZWN0b3I8bG93ZXI9MD5bMl0gc2lnbWFfdGVtcDsKICAgIGNvcnJfbWF0cml4WzJdIFJobzsKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgUmhvIH4gbGtqX2NvcnIoMik7CiAgICBzaWdtYV90ZW1wIH4gZXhwb25lbnRpYWwoMSk7CiAgICAKICAgIGEgfiBub3JtYWwoIDAgLCAuNSApOyAvL25hcnJvdyBwcmlvcnMgdG8gb3ZlcmNvbWUgZGl2ZXJnZW50IHRyYW5zaXRpb25zCiAgICBhcCB+IG5vcm1hbCggMCwgMS41ICk7IAogICAgc2hhcGUgfiBleHBvbmVudGlhbCguNSk7IAp7CiAgICB2ZWN0b3JbMl0gWVlbOF07CiAgICB2ZWN0b3JbMl0gTVU7CiAgICBNVSA9IFthLCBhcF0nOwogICAgZm9yIChqIGluIDE6OCkgWVlbal0gPSBbYV90ZW1wW2pdLCBhcF90ZW1wW2pdXSc7CiAgICBZWSB+IG11bHRpX25vcm1hbChNVSwgcXVhZF9mb3JtX2RpYWcoUmhvLCBzaWdtYV90ZW1wKSk7Cn0KICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBfdGVtcFt0ZW1wc1tpXV07CiAgICAgICAgcFtpXSA9IGludl9sb2dpdChwW2ldKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYV90ZW1wW3RlbXBzW2ldXTsKICAgICAgICAgLy8gYXBwbHkgdGhlIGludmVyc2UgbGluayBmdW5jdGlvbgogICAgICAgIG11W2ldID0gc2hhcGUgKiBleHAoLShtdVtpXSkpOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocFtpXSkgKyBnYW1tYV9scGRmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBnYW1tYSggc2hhcGUsIG11W2ldICk7Cn0KIgptMS43bXYgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuN212LCBkYXRhPWQsIGNoYWlucz00LCBjb3Jlcz00LCBjb250cm9sPWxpc3QoYWRhcHRfZGVsdGE9Ljk5KSkKYGBgIAoKMTAwIGRpdmVyZ2VudCB0cmFuc2l0aW9ucy4uLgoKYGBge3J9CnByZWNpcyhtMS43bXYsIGRlcHRoID0gMykKYGBgCgoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS43bXYpCgpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYV8iCikpICMgdGhlc2UgYXJlIHRoZSBzaGFwZSBhbmQgdGhlIGxvZyhtZWFuKGF2ZyB0aW1lIHRvIGdlcm1pbmF0ZSkpCm11X3JhdGUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBtdXRhdGUoYWNyb3NzKC1zaGFwZSwgfiBzaGFwZSpleHAoLSguKSkpICkgJT4lCiAgc3VtbWFyaXplKGFjcm9zcyhldmVyeXRoaW5nKCksIG1lYW4pKSAlPiUgc2VsZWN0KC1zaGFwZSkgIyBwb3N0ZXJpb3IgbWVhbgptdV9yYXRlCm11X3NoYXBlIDwtIG1lYW4ocG9zdCRzaGFwZSkKCnBvc3RfbG9naXRwIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYXBfIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCgpwb3N0ZXJpb3JfY29lZiA8LSB0aWJibGUoc2hhcGU9bXVfc2hhcGUsIHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIHJhdGU9dChtdV9yYXRlKSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApKQoKcG9zdF9wbG90MS43bXYgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuN212ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCm1ha2UgYSBjb21wYXJpc29uIHBsb3QKCmBgYHtyfQpwb3N0X3Bsb3RfY29tcGFyZSA8LSBiaW5kX3Jvd3Mob2xkPXBvc3RfcGxvdDEuNywgbmV3PXBvc3RfcGxvdDEuN212LCAuaWQgPSAiZml0IikgCnBvc3RfcGxvdF9jb21wYXJlICU+JQogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9c3RyX2ModGVtcHMsZml0KSxsaW5ldHlwZT1maXQpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QsIGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyksIGluaGVyaXQuYWVzID0gRkFMU0UpCiAgCmBgYApvdmVyYWxsIG5vdCBhcyBnb29kCgojIyBleHRlbmQgdG8gbXVsdGlwbGUgc3BlY2llcy9wb3BzCgptb2RpZnkgbW9kZWwgMS43YQoKYGBge3J9Cmxlbmd0aCh1bmlxdWUoZ2VybW9uZSRwb3BzKSkgIzI0IHBvcHMKCmQgPC0gbGlzdChOPW5yb3coZ2VybW9uZSksCiAgICAgICAgICBnZXJtPWdlcm1vbmUkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgcG9wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtb25lJHBvcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMi4xIDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC41ICk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICBzaGFwZSB+IGV4cG9uZW50aWFsKC41KTsgCiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSB7CiAgICAgICAgbXVbaV0gPSBhW3RlbXBzW2ldXTsKICAgICAgICAgLy8gYXBwbHkgdGhlIGludmVyc2UgbGluayBmdW5jdGlvbgogICAgICAgIG11W2ldID0gc2hhcGUgKiBleHAoLShtdVtpXSkpOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocFtpXSkgKyBnYW1tYV9scGRmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBnYW1tYSggc2hhcGUsIG11W2ldICk7Cn0KIgptMS43YSA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS43YSwgZGF0YT1kLCBjaGFpbnM9NCwgY29yZXM9NCwgY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPS45OSkpCmBgYAoKCgoKIyMgaG93IGNhbiB3ZSBzd2l0Y2ggdGhpcyBhcm91bmQgdG8gdXNlIGRheSBhcyBhIHByZWRpY3Rvcj8K